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56230-589 ANAK-244 

DECOMPOSITION OF MULTI-ENERGY SCAN PROJECTIONS 
USING MULTI-STEP FITTING 

Cross References to Related Applications 

[0001 ] There are no prior related patent applications. 

Statement of Government Interest 

[0002] The U.S. Government has no interest in or to the present invention. 

5 Field of the Invention 

[0003] The present invention relates to systems and methods for detecting explosive 

materials using X-ray radiation transmission and scattering to determine one or more physical 
characteristics of a material, and more particularly, to systems and methods for performing 
decomposition of scan projections used in such systems and methods. 

10 Background 

[0004] Various X-ray baggage scanning systems are known for detecting the presence of 

explosives and other prohibited items in baggage, or luggage, prior to loading the baggage onto a 
commercial aircraft. A common technique of measuring a material's density is to expose the 
material to X-rays and to measure the amount of radiation absorbed by the material, the 
1 5 absorption being indicative of the density. Since many explosive materials may be characterized 
by a range of densities differentiable from that of other items typically found in baggage, 
explosives are generally amenable to detection by X-ray equipment. 

[0005] Most X-ray baggage scanning systems in use today are of the "line scanner" type 

and include a stationary X-ray source, a stationary linear detector array, and a conveyor belt for 
20 transporting baggage between the source and detector array as the baggage passes through the 
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scanner. The X-ray source generates an X-ray beam that passes through and is partially 
attenuated by the baggage and is then received by the detector array. During each measuring 
interval the detector array generates data representative of the integral of density of the planar 
segment of the baggage through which the X-ray beam passes, and this data is used to form one 
5 or more raster lines of a two-dimensional image. As the conveyor belt transports the baggage 
past the stationary source and detector array, the scanner generates a two-dimensional image 
representative of the density of the baggage, as viewed by the stationary detector array. The 
density image is typically displayed for analysis by a human operator. 

[0006] Techniques using dual energy X-ray sources are known for providing additional 

10 information about a material's characteristics, beyond solely a density measurement. Techniques 
using dual energy X-ray sources involve measuring the X-ray absorption characteristics of a 
material for two different energy levels of X-rays. Depending upon the calibration of the 
scanner, dual energy measurements provide an indication of dual parameters of the material 
being scanned; for example, at one calibration setting, the dual parameters can be chosen to be 

15 the material's atomic number and the material's density. At another calibration setting, the dual 
parameters can be chosen to be the material's Photoelectric coefficients and the material's 
Compton coefficients. At yet another calibration setting, the dual parameters can be chosen to be 
an amount of a first material present (e.g., plastic) and an amount of a second material present 
(e.g., aluminum). Dual energy X-ray techniques for energy-selective reconstruction of X-ray 

20 Computer Tomography (hereinafter referred to as CT) images are described, for example, in 
Robert E. Alvarez and Albert Macovski, "Energy-selective Reconstructions in X-ray 
Computerized Tomography", Phys. Med. Biol. 1976, Vol. 21, No. 5, 733-744; and United States 
Patent Nos. 4,029,963 and 5,132,998. One algorithm used to generate such dual parameters from 
dual energy X-ray projection data is known as the Alvarez/Macovski Algorithm (hereinafter 

25 referred to as AMA). 

[0007] One proposed use for such dual energy techniques has been in connection with a 

baggage scanner for detecting the presence of explosives in baggage. Explosive materials are 
generally characterized by a known range of atomic numbers and are therefore amenable to 
detection by such dual energy X-ray sources. One such dual energy source is described in U.S. 
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Patent No. 5,661,774, entitled "Improved Dual Energy Power Supply." 
[0008] Plastic explosives present a particular challenge to baggage scanning systems 

because, due to their moldable nature, plastic explosives may be formed into geometric shapes 
that are difficult to detect. Most explosives capable of significantly damaging an aircraft are 
5 sufficiently large in length, width, and height so as to be readily detectable by an X-ray scanner 
system regardless of the explosive's orientation within the baggage. However, a plastic 
explosive powerful enough to damage an aircraft may be formed into a relatively thin sheet that 
is extremely small in one dimension and is relatively large in the other two dimensions. The 
detection of plastic explosives may be difficult because it may be difficult to see the explosive 

10 material in the image, particularly when the material is disposed so that the thin sheet is 
perpendicular to the direction of the X-ray beam as the sheet passes through the system. 
[0009] Thus, detection of suspected baggage requires very attentive operators. The 

requirement for such attentiveness can result in greater operator fatigue, and fatigue as well as 
any distractions can result in a suspected bag passing through the system undetected. 

15 Accordingly, a great deal of effort has been made to design a better baggage scanner. Such 
designs, for example, have been described in U.S. Patent Nos. 4,759,047 (Donges et ah); 
4,884,289 (Glockmann et ah); 5,132,988 (Tsutsui et ah); 5,182,764 (Peschmann et ah); 
5,247,561 (Kotowski); 5,319,547 (Krug et ah); 5,367,552 (Peschmann et ah); 5,490,218 (Krug et 
ah) and German Offenlegungsschrift DE 31 503 06 Al (Heimann GmbH). 

20 [0010] At least one of these designs, described in U.S. Patent Nos. 5,182,764 

(Peschmann et al.) and 5,367,552 (Peschmann et al.) (hereinafter the 764 and '552 patents), has 
been commercially developed and is referred to hereinafter as the "Invision Machine." The 
Invision Machine includes a CT scanner of the third generation type, which typically includes an 
X-ray source and an X-ray detector system secured respectively to diametrically opposite sides 

25 of an annular-shaped platform or disk. The disk is rotatably mounted within a gantry support so 
that in operation the disk continuously rotates about a rotation axis while X-rays pass from the 
source through an object positioned within the opening of the disk to the detector system. 
[0011] The detector system can include a linear array of detectors disposed as a single 

row in the shape of a circular arc having a center of curvature at the focal spot of the X-ray 
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source, i.e., the point within the X-ray source from which the X-rays emanate. The X-ray source 
generates a fan shaped beam, or fan beam, of X-rays that emanates from the focal spot, passes 
through a planar imaging field, and is received by the detectors. The CT scanner includes a 
coordinate system defined by X-, Y- and Z-axes, wherein the axes intersect and are all normal to 
5 one another at the center of rotation of the disk as the disk rotates about the rotation axis. This 
center of rotation is commonly referred to as the "isocenter." The Z-axis is defined by the 
rotation axis and the X- and Y-axes are defined by and lie within the planar imaging field. The 
fan beam is thus defined as the volume of space defined between a point source, i.e., the focal 
spot, and the receiving surfaces of the detectors of the detector array exposed to the X-ray beam. 

10 Because the dimension of the receiving surfaces of the linear array of detectors is relatively small 
in the Z-axis direction the fan beam is designed to be relatively thin in the Z-axis direction. Each 
detector generates an output signal representative of the intensity of the X-rays incident on that 
detector. Since the X-rays are partially attenuated by all the mass in their path, the output signal 
generated by each detector is representative of the density of all the mass disposed in the imaging 

1 5 field between the X-ray source and that detector. 

[0012] As the disk rotates, the detector array is periodically sampled, and for each 

measuring interval each of the detectors in the detector array generates an output signal 
representative of the density of a portion of the object being scanned during that interval. The 
collection of all of the output signals generated by all the detectors in a single row of the detector 

20 array for any measuring interval is referred to as a "projection," or equivalently as a 'View," and 
the angular orientation of the disk (and the corresponding angular orientations of the X-ray 
source and the detector array) during generation of a projection is referred to as the "projection 
angle." At each projection angle, the path of the X-rays from the focal spot to each detector, 
called a "ray," increases in cross section from a point source to the receiving surface area of the 

25 detector, and thus is thought to magnify the density measurement because the receiving surface 
area of the detector area is larger than any cross sectional area of the object through which the 
ray passes. 

[0013] As the disk rotates around the object being scanned, the scanner generates a 

plurality of projections at a corresponding plurality of projection angles. Using well known 
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algorithms a CT image of the object may be generated from all the projection data collected at 
each of the projection angles. The CT image is representative of the density of a two 
dimensional "slice" of the object through which the fan beam has passed during the rotation of 
the disk through the various projection angles. The resolution of the CT image is determined in 
5 part by the width of the receiving surface area of each detector in the plane of the fan beam, the 
width of the detector being defined herein as the dimension measured in the same direction as the 
width of the fan beam, while the length of the detector is defined herein as the dimension 
measured in a direction normal to the fan beam parallel to the rotation or Z-axis of the scanner. 
In general, the resolution of the CT image is inversely proportional to the width of the receiving 

1 0 surface of each detector in the plane of the fan beam. 

[0014] FIG.s 1, 2 and 3 show perspective, end cross-sectional and radial cross-sectional 

views, respectively, of a typical baggage scanning system 100, which includes a conveyor 
system 1 10 for continuously conveying baggage or luggage 1 12 in a direction indicated by arrow 
114 through a central aperture of a CT scanning system 120. The conveyor system includes 

15 motor driven belts for supporting the baggage. Conveyer system 1 10 is illustrated as including a 
plurality of individual conveyor sections 122; however, other forms of conveyor systems may be 
used. 

[0015] The CT scanning system 120 includes an annular shaped rotating platform, or 

disk, 124 disposed within a gantry support 125 for rotation about a rotation axis 127 (shown in 

20 FIG. 3) that is preferably parallel to the direction of travel 1 14 of the baggage 112. Disk 124 is 
driven about rotation axis 127 by any suitable drive mechanism, such as a belt 116 and motor 
drive system 118, or other suitable drive mechanism, such as the one described in U.S. Patent 
No. 5,473,657 issued December 5, 1995 to Gilbert McKenna, entitled "X-ray Tomographic 
Scanning System," which is assigned to the present assignee and which is incorporated herein in 

25 its entirety by reference. Rotating platform 124 defines a central aperture 126 through which 
conveyor system 110 transports the baggage 112. 

[0016] The system 120 includes an X-ray tube 128 and a detector array 130 which are 

disposed on diametrically opposite sides of the platform 124. The detector array 130 can be a 
two-dimensional array such as the array described in U.S. Patent No. 6,091,795 entitled, "Area 
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Detector Array for Computed Tomography Scanning System." The system 120 further includes 
a data acquisition system (DAS) 134 for receiving and processing signals generated by detector 
array 130, and an X-ray tube control system 136 for supplying power to, and otherwise 
controlling the operation of, X-ray tube 128. The system 120 is also preferably provided with a 
5 computerized system (not shown) for processing the output of the data acquisition system 134 
and for generating the necessary signals for operating and controlling the system 120. The 
computerized system can also include a monitor for displaying information including generated 
images. System 120 also includes shields 138, which may be fabricated from lead, for example, 
for preventing radiation from propagating beyond gantry 125. 

10 [0017] The X-ray tube 128 may generate a pyramidically shaped beam, often referred to 

as a "cone beam," 132 of X-rays that pass through a three dimensional imaging field, through 
which conveying system 110 transports baggage 112. After passing through the baggage 
disposed in the imaging field, detector array 130 receives cone beam 132 and generates signals 
representative of the densities of exposed portions of baggage 1 12. The beam therefore defines a 

15 scanning volume of space. Platform 124 rotates about its rotation axis 127, thereby transporting 
X-ray source 128 and detector array 130 in circular trajectories about baggage 112 as the 
conveyor system 110 continuously transports baggage through central aperture 126, so as to 
generate a plurality of projections at a corresponding plurality of projection angles. 
[0018] Pre-reconstruction analysis, post-reconstruction analysis and multiple 

20 projection/non-CT analysis are three prior art techniques generally recognized for using dual 
energy X-ray sources in materials analysis (e.g., in a baggage scanner for detecting the presence 
of explosives in baggage). In pre-reconstruction analysis, the signal flow is as shown in FIG. 4. 
The scanner 120 is typically similar to the one shown in FIG. 1 and has an X-ray source capable 
of producing a fan beam at two distinct energy levels (i.e., dual energy). The DAS 134 gathers 

25 signals generated by detector array 130 at discrete angular positions of the rotating platform 124, 
and passes the signals to the pre-processing element 206. The pre-processing element 206 re- 
sorts the data it receives from the DAS 134 in order to optimize the sequence for the subsequent 
mathematical processing. The pre-processing element 206 also corrects the data from the DAS 
134 for detector temperature, intensity of the primary beam, gain and offset, and other 
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deterministic error factors. Finally, the pre-processing element 206 extracts data corresponding 
to high-energy views and routes it to a high energy channel path 208, and routes the data 
corresponding to low-energy views to a low energy path 210. 

[0019] The projection computer 212 receives the projection data on the high energy path 

5 208 and the low energy path 210 and performs Alvarez/Macovski Algorithm processing to 
produce a first stream of projection data 214 which is dependent on a first parameter of the 
material being scanned and a second stream of projection data 216 which is dependent on a 
second parameter of the material scanned. The first parameter is often the atomic number and 
the second parameter is often material density, although other parameters may be selected. A 

10 first reconstruction computer 218 receives the first stream of projection data 214 and generates a 
CT image from the series of projections corresponding to the first material parameter. A second 
reconstruction computer 220 receives the second stream of projection data 216 and generates a 
CT image from the series projections corresponding to the second material parameter. 
[0020] In post-reconstruction analysis, the signal flow is as shown in FIG. 5. As is 

15 described herein for pre-processing analysis, a pre-processing element 206 receives data from a 
DAS 134, performs several operations upon the data, then routes the data corresponding to high- 
energy views to a high energy path 208 and routes the data corresponding to low-energy views to 
a low energy path 210. A first reconstruction computer 218 receives the projection data from the 
high-energy path 208 and generates a CT image corresponding to the high-energy series of 

20 projections. A second reconstruction computer 220 receives the projection data from the low- 
energy path 210 and generates a CT image corresponding to the low-energy series of projections. 
A projection computer 212 receives the high energy CT data 222 and the low-energy CT data 
224 and performs AMA processing to produce CT data 226 which is dependent on a first 
parameter of the material being scanned and a second stream of projection data 228 which is 

25 dependent on a second parameter of the material scanned. 

[0021] Various approaches have been used for decomposition of projection data (P). For 

example, the AMA non-linear equation approximates the integral equations by a second order 
power series in A c and A py where A c represents the Compton line integral and A p represents the 
Photoelectric line integral. The line integral equations are given as: 
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Pl =b 0 + b x A c + b 2 A p +b z A c A p +b A A] + b s A 2 p (1) 
P H = c 0 + c x A c +c 2 A p + c,A c A p + c A A 2 c + c 5 A 2 p 

[0022] The coefficients b g - and c,- are determined through a calibration procedure as 

5 follows. By measuring the projections values of the combination of various thickness of two 
known materials, a set of equations in the form of Equation 1 can be formed. Since P L , P H , A c 
and A p are known for the calibration data, the coefficients are calculated using a polynomial least 
squares fitting algorithm. 

[0023] Once the coefficients 6, and c, are determined, the decomposition of the Compton 

10 (A c ) and Photoelectric (A p ) images from projections Pl, Ph is accomplished by solving Equation 
1 using the Newton-Raphson method for each pixel in the image. 

[0024] The disadvantage of the method is that it is computationally slow and potentially 

unstable because the iteration can diverge if the initial guess to the solution of Equation 1 is not 
sufficiently close to the true solution. Therefore, the method is susceptible to noise and to 

15 extrapolation outside the region of calibration causing large approximation errors. 

[0025] Another prior art method of performing decomposition is the direct 

approximation method, discussed in L. A. Lehmann, R. E. Alvarez, A. Macovski, W. R. Brody, 
N. J. Pelc, S. J. Riederer, and A. L. Hall, Generalized Image Combinations In Dual KVP Digital 
Radiography, Med. Phys. 8, 659-667 (1981). In the direct approximation method, A c and A p are 

20 expressed as a power series in Pl, Ph 

4=4>+ d x P L + d 2 P u + d 3 P h P H + d 4 P? + d 5 P> (2) 
A P = e o + *A + *2 P h + *APh + erf + 

[0026] As in the non-linear equation method, the coefficients d, and are determined 

25 through a calibration procedure, as follows. By measuring the projections values of the 

combination of various thicknesses of two known materials, one obtains a set of equations in the 

form of Equation 1. Since P L , Ph, A c and A p are known for the calibration data, the coefficients 

are calculated using a polynomial least squares fitting algorithm. For a given pair of projections 
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(Pl, Ph), A c and A p are can be calculated directly from Equation 2. 

[0027] The direct decomposition method avoids the issues of iterative convergence of the 

non-linear equations. However, the method has the disadvantage that it requires a large number 
of polynomial terms in Equation 2 in order to obtain accuracy over an extended range of 
5 projection values. The use of higher order polynomials then introduces a ripple which can be 
controlled only by obtaining a large number of calibration points. 

[0028] In yet another prior art method, decomposition is accomplished using iso- 

transmission lines, described K. Chuang and H. K. Huang, A Fast Dual-Energy Computational 
Method Using Isotransmission Lines and Tables, Med. Phys. 14, 186-192 (1987). According to 
10 this method, for a given projection value, an iso-transmission line is represented by a linear 
equation in two basis functions, wherein the thickness of aluminum and plastic (as examples) is 
given by t AL and tp/, respectively, as: 

PL = at Al + bt Pi (3) 
P H = dt Ai + etpi 

15 where P L and Pu are the high and low energy projections and a, b, d, e are regression coefficients 
that are proportional to the attenuation coefficients of aluminum and plastic. The solution of 
Equation 3 determines the aluminum and plastic combination (t Al , t Pi ) that determines (P L , ^h). 
[0029] The regression coefficients are generated from high and low energy calibration 

tables. The tables contain regression coefficients generated for pre-defined projection values. 

20 The pre-defined projections are generated by scanning known thicknesses of aluminum and 
plastic. For a given projection value, the regression coefficients are generated by interpolation 
between the closest two pre-defined coefficients stored in the table. 

[0030] The iso-transmission line method requires a large amount of calibration data in 

order to generate the regression coefficients a, b, d, e. Further, the iso-transmission lines become 
25 increasingly non-linear as the projection value increases. In such a situation, the linear equations 
are not valid and the method causes large approximation errors. 

[0031] Other prior methods for dual energy decomposition also suffer from 

shortcomings, as discussed in H. N. Cardinal and A, Fenster, An Accurate Method For Direct 
Dual Energy Calibration and Decomposition by, Med. Phys. 17, 327-341 (1990). As an 
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example, the direct decomposition algorithm produces incomplete separation of objects, as 
determined visually. Additionally, the direct decomposition algorithm also introduces artifacts 
in the decomposed images of scanner data. The artifacts are primarily produced by 
approximation errors, sensitivity to noise and the lack of handling of exceptions in the input 
5 projection data. In addition, the algorithm requires a calibration procedure that is complicated, 
time consuming and prone to measurement errors, which propagate through the decomposition 
procedure to be manifested as artifacts in images. 

Summary of the Invention 

10 [0032] In accordance with the present invention, a dual energy decomposition algorithm 

with multi-step fitting is provided that better separates objects with fewer artifacts and produces 
fewer false alarms. In accordance with the present invention the dual energy decomposition 
algorithm may be comprised of two sub-processes or parts (or sub-algorithms). First, the 
decomposition algorithm that operates on received or detected scanner data. Second, if not yet 

15 performed, prior to performing the decomposition, a calibration procedure may be performed to 
calibrate the decomposition algorithm. Accordingly, below, the dual energy decomposition 
algorithm is discussed and the optional calibration procedure is also discussed. Generally 
speaking, these two sub-process execute independently, at different times, so are discussed 
separately. 

20 [0033] The calibration procedure can use simulated data or measured data or a 

combination of simulated and measured data. Using simulated data or stored measured data 
allows the calibration procedure to be relatively quick to perform and free from noise. 
Preferably, such a calibration procedure is not machine-specific and can be done once per 
scanner design. In still other embodiments, a noise reduction module or step may be included 

25 prior to processing projection data using the dual energy decomposition algorithm. 

[0034] Dual energy decomposition with multi-step fitting in accordance with the present 

invention is useful for the decomposition of projection data, wherein such projection data may 
include input projection data acquired using at least two x-ray spectra for a scanned object, 
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including low energy projection data (Pi) and high energy projection data (Ph). The method 
comprises solving the projections Pi and P\\ to determine a photoelectric line integral (A p ) 
component of attenuation and a Compton line integral (A c ) component of attenuation of the 
scanned object using a multi-step fitting procedure and reconstructing a Compton image I c and a 
5 photoelectric image I p from the Compton line integral and photoelectric line integral. 

[0035] The method may further include, prior to solving the projections, performing a 

calibration procedure using simulated data, wherein the calibration procedure may include 
generating low energy iso-transmission contours for known values of Pi at known values of the 
photoelectric line integral A p and at known values of the Compton line integral A c . According to 

10 the method, for each low-energy iso-transmission contour corresponding to P L , minimum and 
maximum values of P H may be computed, i.e., PH_min and P H _max- A polynomial m L may be used 
to determine a minimum bound on P H > and the multi-step fitting procedure may further include a 
fitting procedure for fitting the minimum values of P H to a polynomial function mi(P L ). Also, 
the method may further include using a polynomial n L to determine a maximum bound on P H , 

15 and the multi-step fitting procedure may include a fitting procedure for fitting the maximum 
values of Ph to a polynomial function ni(P\). 

[0036] The multi-step fitting procedure may include a fitting procedure wherein, for each 

of the low energy iso-transmission contours, A p is fit to a polynomial function gi(A c ). And, the 
calibration procedure may further include fitting the set of coefficients gu determined for known 

20 values of P L to a polynomial function hu(Pi). 

[0037] The calibration procedure may also include generating high energy iso- 

transmission contours for known values of Ph at known values of A p and at known values of A c . 
The multi-step fitting procedure may then include for each of the high energy iso-transmission 
contours, fitting A p to a polynomial function of gn(A c ). And, the calibration procedure may 

25 further include fitting the set of coefficients g H i determined for known values of P H to a 
polynomial function hw\(P\\), 

[0038] Apart from calibration, the dual energy decomposition procedure with multi-step 

fitting method may also include generating a low energy iso-transmission contour corresponding 
to Pi and a high energy iso-transmission contour corresponding to Ph- The method may then 
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include determining the values of the photoelectric line integral A p and Compton line integral A c 
at the intersection of the low energy iso-transmission contour and the high energy iso- 
transmission contour, wherein the intersection of the low energy iso-transmission contour and 
the high energy iso-transmission contour may be determined by equating a first polynomial 
5 function representing the low energy iso-transmission contour with a second polynomial function 
representing the high energy iso-transmission contour. In such a case, the coefficients of the first 
polynomial function may be determined using P L and the polynomial function h L (P L ) generated 
during calibration, and wherein the coefficients of the second polynomial function may be 
determined using P H and the polynomial function h\\ (P H ) generated during calibration. 

10 [0039] The method may include computing a modified value of the input low energy 

projection data and a modified value of the input high energy projection data, wherein each of 
the modified values, designated as Plc for low energy projection data and P Hc for high energy 
projection data, is clamped to be bounded between two values. The coefficients of the low 
energy iso-transmission contour may be computed for the value of Pu using the polynomial 

15 function hi generated during calibration, as described above, wherein hi describes the variation 
of the coefficients of the polynomial gu which describes the shape of the low energy iso- 
transmission contour for the value of P Lc . The coefficients of the high energy iso-transmission 
contour may be computed for the value of P Hc using the polynomial function h H generated during 
calibration, as described above, wherein h H describes the variation of the coefficients of the 

20 polynomial which describes the shape of the high energy iso-transmission contour for the 
value of Pho The modified value P^ is computed by clamping the value of P L to lie between 0 
and the maximum value of Pi used to generate the calibration data. In accordance with the 
method, the modified value Ph c niay be determined by clamping the value of Ph to lie between a 
minimum value of P H (i.e., P H min ) and a maximum value of P H (i.e, P H max ). ^H min and P H ™ X are 

25 determined using the polynomials m L and m f determined during calibration, and as a function of 
Plc, wherein m L is a polynomial used to determine Ph 1 ™ and n L is a polynomial used to 
determine Ph™*. 

[0040] The method may also include calculating a scaled Compton line integral value 

(A cs ) as a function of a scale factor s c and A c and calculating a scaled photoelectric line integral 
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value (A ps ) as a function of a scale factor s p and A ps . Scaling A c and A p is performed to obtain 
pixel intensities within a specific range of values if A c and A p are reconstructed into images. The 
method may further include reconstructing image I c and image I p as a function of A cs and A ps . 
The method may then include determining an image of a basis function X(I X ) and said basis 
5 function Y(/y), by solving image I c and image I p on a pixel-by-pixel basis. The basis functions 
are two functions that can be combined to determine the pixel intensities in the Compton image 
and a photoelectric image I p . 

[0041] A system for decomposing projection data using multi-step fitting that includes 

modules configured to implement the above functionality may also be provided. Such a system 

10 may include, for projection data of a set of scanned objects acquired using at least two x-ray 
spectra, media for storing low energy projection data (P L ) and high energy projection data (P H ). 
The system may include a calibration module configured to calibrate the decomposition module 
using at least some simulated data; generate a low energy iso-transmission contour 
corresponding to known values of Pi and a high energy iso-transmission contour corresponding 

15 to known values of Ph; generate a polynomial gL that represents the shape of the low energy iso- 
transmission contour; and generate a polynomial gu that represents the shape of the high energy 
iso-transmission contour; generate polynomials hi that represent the variation of the coefficients 
of the polynomial gL as a function of Pi\ generating polynomial h H that represents the variation 
of the coefficients of the polynomial g H as a function of P H ; and determining the minimum and 

20 maximum values of Ph for each transmission line corresponding each P L ; and generating a 
polynomial tn\\ that represents the variation of the minimum value of Ph as a function of Pi\ and 
generating a polynomial rtu that represents the variation of the maximum value of Ph as a 
function of P L ; 

[0042] A decomposition module may be provided that is configured to determine a 

25 photoelectric line integral (A p ) component of attenuation and a Compton line integral (A c ) 
component of attenuation for measured values of P L and Ph using multi-step fitting, and 
configured to determine coefficients of gL using P L and hi(Pi) and determine coefficients of gn 
using Ph and hniPn) wherein A c and A p are determined as a function of Pi and P H using the 
coefficients of g L and the coefficients of g H . 
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[0043] The system may also include an image reconstruction module configured to 

reconstruct a Compton image (7 C ) and a photoelectric image (/ p ) from the A p and A c . The image 
reconstruction module may be configured to determine an image of a basis function X(/ x ) and of 
a basis function Y(/ Y ), by solving image I c and image I p on a pixel-by-pixel basis. 

Brief Description of the Drawings 

[0044] The drawing figures depict preferred embodiments by way of example, not by 

way of limitations. In the figures, like reference numerals refer to the same or similar elements. 

[0045] FIG. 1 is a perspective view of a baggage scanning system, known in the prior art. 

[0046] FIG. 2 is a cross-sectional end view of the system of FIG. 1. 

[0047] FIG. 3 is a cross-sectional radial view of the system of FIG. 1. 

[0048] FIG. 4 is a signal flow diagram of a system capable of performing 

pre-construction analysis, useful in the system of FIG. 1. 

[0049] FIG. 5 is a signal flow diagram of a system, like that of FIG. 1, capable of 

performing post-reconstruction analysis. 

[0050] FIG. 6 is a flowchart showing a method of dual energy decomposition with multi- 

step fitting, in accordance with the present invention. 

[0051] FIG. 7 is a vector representation of photoelectric and Compton coefficients for 

aluminum and plastic, in accordance with the illustrative embodiment. 

Detailed Description of the Preferred Embodiments 

[0052] In accordance with the present invention, a dual energy decomposition algorithm 

is provided that better separates objects with fewer artifacts and produces fewer false alarms, and 
preferably uses at least a two step fitting process, generally referred to as "dual energy 
decomposition with multi-step fitting". A decomposition algorithm in accordance with the 
present invention may also use a calibration procedure that is based at least in part, if not 
entirely, on simulated data, and therefore is quick to perform and free from noise often 
introduced by other forms of calibration data. The calibration procedure may also be performed 
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using measured data or a combination of simulated data and measured data. Preferably, such a 
calibration procedure is not machine-specific and can be done once per scanner design. In still 
other embodiments, a noise reduction module or step may be included prior to processing of 
projection data. Given that the calibration procedure and the dual energy decomposition 
algorithm operate relatively independently, they are generally discussed below separately. 
[0053] For the range of photon energies used for medical and commercial diagnostic 

imaging (i.e., from about 20 to 200 keV), the attenuation of an object is determined by the 
photoelectric absorption and Compton scatter. Photoelectric absorption consists of an x-ray 
photon imparting all of its energy to a tightly bound inner electron in an atom. The electron uses 
some of its acquired energy to overcome the binding energy within its shell, the rest appearing as 
kinetic energy of the thus freed electron. Compton scatter consists of the interaction of the x-ray 
photon with either a free electron, or one that is loosely bound in one of the outer shells of an 
atom. As a result of Compton scatter, the x-ray photon is deflected from its original direction of 
travel with some loss of energy, which is gained by the electron. Both the photoelectric and 
Compton effects are energy dependent. 

[0054] The linear attenuation coefficient, fi{\,yJE) for a given material as a function of 

energy E and spatial position (x,y) in a plane can be approximated as: 

M(x,y,E) = a p {x,y) %(E) + a c (x,y) O c (£) (4) 
where <& P (E) is the photoelectric energy dependence given by: 

<M*) = ^T (5) 
and O c (£) is the Compton energy dependence that is described by the Klein-Nishina function 
given by: 

^ )= I±420^_i ln(l+2 J + i ln(l+2 „)_l±M 

a |_(l + 2tf) a '\ 2a v ' (i + 2a) 2 V ' 

where a =E/5 10.975 keV. The parameters a p (x,y) and a c (x,y) are constants that depend on the 
material composition, namely: 

a P (x,y) =p e (x,y)BZ'(x,y) (7) 
ac{x,y)=p e (x,y) 
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where B = 9.8xl0~ 24 is a constant, n * 3, Z(x,y) is the atomic number, and p e (x,y) is the electron 
density given by: 



p e (x,y) = N A 



( Z(x,y) 



(8) 



where N A is Avogadro's number (6.023 x 10 23 ) and W(x,y) is the atomic weight. 
5 [0055] In accordance with the present invention, the Compton coefficient, a c (x,y), and the 

photoelectric coefficient, a p (x,y), of an unknown material may be determined on a per-pixel basis 
from images reconstructed from projection measurements using at least two source spectra, i.e., 
low energy and high energy. A dual energy decomposition algorithm with multi-step fitting, in 
accordance with the present invention, comprises the following steps, described in greater detail 
10 below: 

1) The projections of the objects are measured using at least two different source 
spectra, obtained by varying voltage to the x-ray tube and optionally changing the beam 
filtration. 

2) The projections (P) are decomposed into line integrals of the photoelectric and 
1 5 Compton components of attenuation using a two-step fitting procedure. 

3) The photoelectric and Compton line integrals are filtered and "back projected" to 
reconstruct photoelectric and Compton images. 

4) Optionally, the photoelectric and Compton images are solved on a pixel-by-pixel 
basis to generate images of two known basis functions, such as aluminum and plastic, as one 

20 example, or mass and atomic number, as another example. 

[0056] FIG. 6 is a flow chart 600 showing the above steps of dual energy decomposition 

in accordance with the present invention. Input projections P L 610 (low energy) and P H 620 
(high energy) are decomposed into Compton and photoelectric line integrals, A c 630 and A p 640, 
respectively. A c and A p are back-projected to generate the Compton and Photoelectric images I c 

IS 650 and I p 660. I c and I p are optionally solved to determine the images of the basis functions, 
e.g., for plastic and aluminum, I Pi 670 and I AI 680. "Solving" the images is a mathematical 
operation described in more detail below. 
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I: Overview of Decomposition with Multi-Step Fitting 

[0057] An overview of the mathematics of present dual energy decomposition with 

multi-step fitting is provided below. This embodiment, which is also applicable to three- 
dimensional data, is described with respect to a CT scanner with an x-ray source and detector 
configuration, such as that shown and described with respect to FIG.s 1, 2, and 3. For this 
example, an object with attenuation n{x,y,E) is scanned. An example of one such scanner is the 
AN6000 EDS scanner by Analogic Corporation, Peabody, Massachusetts. The scanner source 
generates a low energy spectrum S L and a high energy spectrum Sr. 4 and I H are the 
corresponding x-ray intensities that are detected when the object is scanned, given by: 



(») 



I oL and 7oh are the corresponding x-ray intensities that are detected in the absence of any object, 
given by: 



l oL 



= [S L (E)dE 



(10) 



I oH = [S L (E)dE 

and lfi(x,y,E)ds is the path integral, which represents the attenuation of x-rays as they travel 
along a straight line through the object. 

[0058] P L and P H are the corresponding low and high energy projections that are 

measured by the scanner's detector, given by: 



P L =-ln 



P H =-ln 



_f_H 

v^oHy 



= -lnf S L (Ey^ x ' y ' E}b dE + ln \S L {E)dE 
= -ln JT S H (Eyl Ax ' y - Eyis dE+\n \S H (E)dE 



(11) 



[0059] 



From Equation 4 for fi{x,y,E), the line integral j^(x,y,E)ds can be written as: 

J ^x,y,E)ds = A p O p (E)+ A c O c {e) (12) 
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where A c and A p are the line integrals of a c (x,y) and a p (x,y) given by: 

Ac = J (13) 

[0060] Substituting Equation 12 into Equation 1 1 yields: 

P, =A{A C ,A P )= -In jS L {E)e- A >«> {EyA <*' {E) dE + ln js t (E)dE (14) 

=/„U,^)=-ln jS H (E)e- A ^ E) - A ^ E) dE + \n \S H (E)dE 

[0061] Equation 14 can then be solved to determine the values (A c , A p ) that generate the 

projections (P L , P H ) provided that the Jacobian (J), given by: 



J=det 



dA c dA p 
dP H dP H 



dA c dA„ . 

C p J 



(15) 



is non-zero, where det represents the determinant. 

[0062] Note that the direct analytical solution of Equation 14 is difficult to obtain 

because the exact analytic forms of the spectra, S L and 5 H , are not easily measured. Therefore, 
several methods are described in the prior art to determine the solution of Equation 14 without 
knowledge of the spectra such as the direct decomposition method or the iso-transmission 
method, previously discussed. According to the present algorithm for dual energy decomposition 
with multi-step fitting, the solution of Equation 14 is preferably determined using the iso- 
transmission method. 

[0063] Accordingly, the iso-transmission method used in the preferred embodiment, a 

given projection value P L can be determined by multiple values of (A c , A p ). The values lie on a 
contour, which is defined as an iso-transmission contour. Each contour is described by the 
equation: 

4>=Sl(4)U (16) 
where g L is a polynomial function that describes the shape of the contour for a given P L given 
by: 
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& (Ac) \ ^=gu> + giJc + guA 2 c + . 



(17) 



10 



where gu is the polynomial coefficient of order i for a given value of P L . 
[0064] Note that the polynomial must satisfy the following properties of the contour. 

The values of A c and A p are greater than or equal to zero and, therefore, lie in the first quadrant of 
a two dimensional coordinate system. The value of A p decreases as the value of A c increases. 
Therefore, for a given A c , there is only one value of A p that determines P L . Therefore, g L is a 
monotonically decreasing function within the first quadrant and is bounded the values (A max, 0) 
and (0, ^p max) where A c mix is the solution of the equation: 

& (^c_max) I ^ = 0 (18) 

and A pmM is determined as: 

4u»x=gL(0)| Pl (19) 



1 5 [0065] The value of P H varies continuously along the contour according to the equation: 

p n\ p L =MA c ,g L (A c )\ Pl ) (20) 

where / H is given by Equation 14. Since g L is a monotonically decreasing function within the 
first quadrant, P H must lie between the following maximum and minimum values, given by: 
20 C'-I^/hKh^O) (21) 

[0066] For a given P H , the iso-transmission contour is described by the equation: 

Ap=gH(A c )\ PH (22) 
where g H is a polynomial function that describes the shape of the contour for a given P H , given 
25 bv: SH(A c )\ PH = g H0 +g Hi A c +g m A 2 c +.... (23) 

where gm is the polynomial coefficient of order j for a given value of P H . 
[0067] In the preferred form, a calibration procedure is used to determine the coefficients 
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of gi and g H for arbitrary projection values. The calibration steps used for the determination of 
gi are as follows. An iso-transmission contour is generated for a known value of Pi using 
known values of (A c , A p ) and simulated high and low energy spectra. The function gi is fit to the 
iso-transmission contour to determine the coefficients of gi for the value of Pi. The coefficients 
5 gu are generated using the above procedure for several known values of Pi up to a maximum 
value of Pijmx- The variation of each coefficient with projection value is fit to a polynomial 
function hi given by: 

gu = *hja + KA + Kh P l + •••• (24) 
where hu\ is the polynomial coefficient of order j (0, 1, 2...) that determine the coefficient gu. 
10 The function h L is used to compute the coefficients of gu for any value of P L . The above 
procedure is repeated to generate the coefficients of the polynomial function ft H > which is used to 
compute the coefficients of gu for any value of P H . 

[0068] Additionally, the variation of the minimum and maximum values of Ph for each 

iso-transmission contour are fit to polynomial functions given by: 
15 P™ \ p =m L0 +m u P h +m L2 Pt +... (25) 

Pr \p L = n L0 + n Ll P L +n u P?+... 

[0069] Once the polynomials gi and g H are determined for the projections (P L , P H ), the 

solution of Equation 14 are determined by the coordinates of the point of intersection of the two 
contours determined by Equations 17 and 22 as: 

20 gL(^c)|p L =gH(^c)|p H (26) 

provided that the values (P L > Ph) satisfy the conditions: 

0<P L <Pr X (27) 

pr \p L <Pn<pr\ PL 

where P™ x is the maximum value of P L used in calibration, and P H mi11 | ^ and P™* | Pi are 

25 determined using Equation 25. 

[0070] Once the values of (A c , A p ) are determined for all projections, (A c , A p ) are back- 

projected to reconstruct a Compton image I c and a photoelectric image 7 P . The intensity of a 
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given pixel (x, y) in I c is the value a c (x, y). The intensity of a given pixel (x, y) in I p is the value 

y)- 

[0071] As an optional step, the photoelectric and Compton images are then solved on a 

pixel-by-pixel basis to generate images of two known basis functions. Examples of basis 
5 functions are the amount of aluminum and plastic, or mass and atomic number. Continuing with 
the illustrative embodiment, the basis functions for aluminum and plastic are used. One reason 
for using aluminum and plastic for illustrative purposes is that the Compton and photoelectric 
coefficients of aluminum and plastic differ from each other significantly (see Table 1) and span a 
wide range of atomic numbers that are encountered in baggage scanning. 

10 



Coefficient 


Aluminum 


Plastic 


Units 




0.390 


0.1952 


cm' 1 


a p 


69734.0 


3309.0 


keV 3 /cm 



Table 1: Compton & photoelectric coefficients for aluminum and plastic 



[0072] The aluminum and plastic images are generated based on the following theory, 

15 demonstrated by the vector representation 700 shown in FIG. 7. A given pixel location (e.g., x, 
y) in I c and I p determines the coefficients a c (x, y), a p (x, y), which can be represented in polar 
coordinates (r 712, 0 714), as the vector V 710 in the two-dimensional space of Compton 
coefficients, a C9 and photoelectric coefficients, a p9 where: 

r(V) = ^a c {x,y) 2 +a p {x,y) 2 (28) 
20 ^(V) = tan- 1 f^^ > 

{"My)) 

[0073] In Equation 28, aluminum is represented by the vector V A i 720 represented in 

polar coordinates as (r A1 722, 0 M 724) and plastic is represented by the vector V P i 730 represented 
in polar coordinates as (r P i 734, 0 Pi 734) in FIG. 7. V A i and V P i are defined to be the unit vectors 
along the directions of V A i and V P] . The image pair (7 C , I p ) can be processed into a pair of 
25 aluminum and plastic images, (7 C , I p ) wherein vectors determined by pixels in (7 C , I p ) are 
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represented as admixtures of aluminum and plastic given by: 



V = 1 V +/ V 



(29) 



where 



I Al = rsin 





(0- 






f 
















2 









(30) 



and #ai and #pi are determined from the Compton and photoelectric coefficients of aluminum and 
plastic using Equation 28. 



II. Details of the Decomposition 

[0074] This section provides specific mathematical details of a dual energy 

10 decomposition algorithm with multi-step fitting, in accordance with the present invention. The 
inputs, parameters and outputs of the algorithm are listed, and then the mathematical details of 
the algorithm are presented below. 

[0075] The inputs of the decomposition algorithm with multi-step fitting include low 

energy projection Pi and high energy projection Pu, and the outputs include aluminum image /ai 
15 and plastic image /pi, in the illustrative embodiment. The various parameters are as follows: 



Symbol 


Definition 




Fit coefficients used to determine low energy iso-transmission contours 




Fit coefficients used to determine high energy iso-transmission contours 


"*Lj 


Fit coefficients used to determine the minimum value of Pr for a given P L 


20 n Lj 


Fit coefficients used to determine the maximum value of Ph for a given Pi 


„ pi 


Compton coefficient of Plastic 


„ pl 

0p 


Photoelectric coefficient of Plastic 


„ Al 


Compton coefficient of Aluminum 


„ Al 


Photoelectric coefficient of Aluminum 


25 s c 


scale for decomposed Compton line integrals 
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Sp scale for decomposed photoelectric line integrals 

k ap Calibration scale factor for A p 



10 



[0076] The preferred decomposition algorithm includes the following steps: 

[0077] Step 1: Where P L is the input low energy projection data, compute the clamped 

value Plc as: 

' 0 P L (0 

p u =i p L o<p L <pr 01) 
ft p L )pr 

[0078] Step 2: Given the value of P u , calculate the values P H min and / > H max as: 



pmin 
H 

nmax 



(32) 



[0079] Step 3: Where P H is the input high energy projection data, compute the clamped 

value P Hc as: 



Phc = 



nmax 



p*)pr(Pu) 



(33) 



[0080] Step 4: Given the input coefficients huj, generate the coefficients gu of the low 

1 5 energy iso-transmission contour as: 



7=6 



y=0 



(34) 



[0081] Step 5: Given the input coefficients /*Hij, generate the coefficients gu\ of the high 

energy iso-transmission contour as: 



(35) 



20 [0082] Step 6: Determine the intersection of the iso-transmission contours corresponding 

to Pic and Ph c as the solution to the equations: 
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as follows: 



£L0 + SLX A c + SllA = £h0 + 8hA + ( 36 ) 

1 ) Calculate the discriminant D as: 

5 2) Clamp the value of D to A calculated as: 

[D D>0 

D=\ (38) 

c [0 D< 0 

This step is performed because, in some cases, negative values of D could be 
calculated due to noise and numerical errors in the data. 

3) Determine the value of A c at the intersection of the high and low energy 
10 contours as the smaller of the roots determined by the equation: 

Ac= -(gu-g w )±f; (39) 

2 UL2-SH2) 

Note that the negative solution from the above equation is clamped to 0 because 
A c is a line integral, which must be non-negative. 

4) Calculate the value of A p using the equation: 

15 A p =g L0 +g u A c + gL2 A 2 c (40) 

5) Calculate A pc as the clamped value of A p : 

^ BC =\ / n (41) 

[0083] Step 7: Given the scale factors s c and s p , calculate the scaled projections as A cs 

and A ps : 

A — s~A 

20 c (42) 

dps - S p^pc 

[0084] Step 8: Reconstruct Compton and photoelectric images I c and I p , from the values 

A cs and A ps using filtered back-projection. 

[0085] Step 9: This is an optional step. Solve the images I c and I p on a pixel-by-pixel 
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basis to determine 7 A i and I ?h as previously described. The equations are repeated here for 
convenience. Let (x, y) be the location of a pixel. The intensity of the pixel in the images (/ A1 , 
/pi) are calculated as: 

r 0(x,y)-0 p ^ 



I P ,{x,y) = r(x,y)cos 
I Ax,y)=r(x,y)sm 



r 0{x,y)-e p ; 



(43) 



where 



(44) 



and 



9 n = tan" 



0 M = tm~ l 



* P op 
V S c a c ) 

r P dp 
\ c c 



(45) 



III. Overview of Calibration 

[0086] This section provides the details of a calibration procedure, which generates 

parameters use by the dual energy decomposition algorithm with multi-step fitting. The 
calibration procedure includes the following steps. 

[0087] Step 1 : Simulated projection data are generated for known values of high and low 

energy projections. 

[0088] Step 2: An iso-transmission contour is determined for each projection. 

[0089] Step 3: For each iso-transmission contour, the photoelectric line integral is fit 

using a quadratic function of the Compton line integral. 

[0090] Step 4: Fit coefficients from the quadratic function are determined for several 

projection values. 

[0091] Step 5: The variation of each fit coefficient with projection value is described by 



25 



BST99 1355645-3.056230.0589 



a polynomial function. 

[0092] Step 6: The fit coefficients of the polynomial function are stored for high and low 

energy values and input to the decomposition algorithm, in accordance with the present 
invention. 

5 [0093] The inputs of the calibration procedure are the simulated low energy input 

spectrum Si and the simulated high energy input spectrum Sh. The outputs include the fit 
coefficients huj used to determine the low energy iso-transmission contours, the fit coefficients 
A H ij used to determine the high energy iso-transmission contours, the fit coefficients w L j used to 
determine the minimum value of Pu for a given and the fit coefficients « L j used to determine 
10 the maximum value of P H for a given P L . 

[0094] The parameters of the calibration procedure are as follows: 

Symbol Definition 

flj 1 Compton coefficient of plastic 

a ? p l Photoelectric coefficient of plastic 

1 5 Nj S0 Number of points generated per iso-transmission line 

Minimum thickness of plastic used in calibration 
t™* Maximum thickness of plastic used in calibration 

dp\ Thickness increment 

k ap Calibration scale factor for A p 

20 

[0095] Below is described the mathematics for generating iso-transmission contours for 

low-energy projections. The same procedure is used to generate iso-transmission contours for 
high-energy projections. 

[0096] Step 1: Given the parameters /™ n , /™ x and d ?h compute the number of projection 

25 values follows: 

.max _ ,min 

"pi 

[0097] Step 2: For each value of t n [n]: 
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'p.M=C + «4>.» » = 0,....,» pn) j (47) 
compute the low-energy projection P L [n] as follows, using Equation 14: 

PM=A(aX[n\a%[n]) (48) 
[0098] Step 3: For each value of the low-energy projection P L [n], compute the iso- 

transmission contour. The iso-transmission contour is a set of ordered pairs of (A c [l][n], 
where / = 0, . .. , N iso , all of which produce the value P L [n]. There are generally three 
steps involved in computing each iso-transmission contour: 

[0099] 1) Compute ^ c [0][«], which is the maximum value of A c for low-energy 

projection P L [«], by iteratively solving the following equation: 

fM0][40) = P L [n] (49) 
[001 00] 2) Compute A c [l\ [n] for each / as follows: 

X[o]M) 



A MM = 4 MM- 1 



^r^' /=i, (so) 

V iy iso J 

[00101] 3) Compute the associated A p [l\[n] by iteratively solving the following 
equation: 

AkMM^MMM /=u..,tf to (si) 

Note that A p [0][n]=0. Below is a more detailed description of the preferred iterative method used 
to solve Equations 49 and 5 1 . 

[00102] 4) Each A p [l\[n] is scaled as follows: 

^J*]=4,[*R, (52) 



where k ap is the scale factor provided to the calibration procedure. The scale factor k ap is used to 
obtain similar ranges of values among A p and A c . 

[00103] Step 4: For each iso-transmission contour (given n), there are N t ordered pairs of 
(A c [[\[n), A p [t\[n\). The N, ordered pairs are fit using the quadratic function: 

A pk = §L0 M+ Su Wc + Su Wl (53) 

[00104] Step 5: Fit the coefficients of the quadratic function using sixth-order polynomial 
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functions of the low-energy projection value P L : 

6 

Su=XVl« i ' = 0 > 1 »2 (54) 

y=o 

The fitting is performed using a chi-squared minimization routine. This routine is known in the 
mathematical arts, and is described in W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. 
Flannery, Numerical Recipes in C, The Art of Scientific Computing, 2nd Ed. (Cambridge 
University Press, New York, NY, 1992), Chap 15, pp. 681-706. 

[00105] Step 6: In addition to the iso-transmission contour, compute the minimum and 
maximum values of the corresponding high-energy projection P^ n [n] and P™*[n] for each 
Pi[n], as follows: 

= /„(4[o],o) 

PTM = /„(<U,M (55) 

[00106] Step 7: Fit the minimum and the maximum P H to sixth-order polynomial 
functions of P L as follows: 



Pt^rn^i 



(56) 



7-0 

[00107] The steps (1) through (6) above is repeated for obtaining coefficients h m for 
generating the iso-transmission contours of the high-energy projections P H . 
[00108] Step 8: The fit coefficients h U j, h m> w L j, and « L j are stored for use by the 
decomposition algorithm. 

[00109] The solution of Equations 49 and 51 is determined using the following iterative 
algorithm. The algorithm is described for a single low energy iso-transmission curve or contour 
for projection P 0 , for a single value of the Compton line integral A c . The same procedure is 
applied to determine values for the high-energy iso-transmission curve or contour. 
[001 10] Step 1: For the value of the Compton line integral A c , two initial estimations for 
the photoelectric line integral are made to obtain projections with magnitudes greater than and 
less than P 0 . Let the estimations be designated A p \ and A p2 , respectively, wherein: 
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A — A. 

A P 2 = "5 

Note that the negative value of A p i is chosen to obtain a projection value that is greater than P 0 . 
The positive value of A p \ = A c is chosen to obtain a projection value that is lesser than P 0 . 
[001 1 1] Step 2: The following steps are iterated. 
5 [001 12] 1) For iteration i, the values of the corresponding low-energy projections are 
calculated using Equation 14. 

pL=f l {a c ,aJ 

where the superscript i is used to represent values from the i iteration and: 

A-* H" A 

^=- £L y- ei (59) 
10 [00113] 2) For iteration /, calculate new values for A' pl and A' p2 as follows: 

ifPL' l <P 0 

pm J m o 




otherwise ^ ^ 

pl otherwise 

[001 14] Step 3: Repeat the Steps 2(1) and 2(2) until the following condition is satisfied: 

Pl-P}< e (62) 
1 5 where e = 0. 1 is the tolerance for convergence. 
[00 1 1 5] The value of A l pm is then returned. 

[001 16] While the foregoing has described what are considered to be the best mode and/or 
other preferred embodiments, it is understood that various modifications may be made therein 
and that the invention or inventions may be implemented in various forms and embodiments, and 
20 that they may be applied in numerous applications, only some of which have been described 
herein. For example, reference was made variously to use of polynomial equations, quadratic 
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equations and equations having variables raised to other powers. However, in other 
embodiments, appropriate equations could take other forms, such as any of a number of forms of 
parametric equations. That is, as will be appreciated by those skilled in the mathematical arts, 
there can be various manners of representing or modeling physical, electrical, mechanical and 
5 other sciences and the equations used to do so can be made more or less complex given whatever 
assumptions are made with respect to the parameters related to such equations. 
[00117] As used herein, the terms "includes" and "including" mean without limitation. It 
is intended by the following claims to claim any and all modifications and variations that fall 
within the true scope of the inventive concepts. 
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